/*-
 * Copyright (c) 2013 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ MTBLAS TEST: tst_lv_3e_gmm.c    Rev. 0.0.4    2013-11-15T16:54+08:00 $
 */

#include <stdlib.h>
#include <stdio.h>
#include <complex.h>

#include <time.h>

#include <cblas.h>

#include <tst_lv_a_com.h>

static const int flg_trn[3] = {MTBLAS_TRANS_NONE, MTBLAS_TRANS, MTBLAS_TRANS_CONJ};


int stblas_3e_mxmxmulacc_test(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double *ma, const size_t lda, const double *mb, const size_t ldb, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi, zj, zk;

    double buf_ma;

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulacc_test()\n");
#endif

    for (zi = 0; zi < rz; zi++)
    {
        for (zk = 0; zk < kz; zk++)
        {
            buf_ma = ma[zi * lda + zk];

            for (zj = 0; zj < cz; zj++)
            {
                mc[zi * cz + zj] += buf_ma * mb[zk * ldb + zj];
            }
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulacc_test() PASSED\n");
#endif

    return ierr;
}


int stblas_3e_mxmxmulacc_test_v1(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double *ma, const size_t lda, const double *mb, const size_t ldb, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;

    double buf_ma;
    double *ptr_ma_0, *ptr_ma_1;
    double *ptr_mb_0, *ptr_mb_1;
    double *ptr_mc_0, *ptr_mc_1;

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulacc_test_v1()\n");
#endif

    ptr_mc_0 = mc;

    for (ptr_ma_0 = ma; ptr_ma_0 < ma + rz * lda; ptr_ma_0 += lda, ptr_mc_0 += ldc)
    {
        ptr_ma_1 = ptr_ma_0;

        for (ptr_mb_0 = mb; ptr_mb_0 < mb + kz * ldb; ptr_mb_0 += ldb)
        {
            buf_ma = *ptr_ma_1;
            ptr_mb_1 = ptr_mb_0;

            for (ptr_mc_1 = ptr_mc_0; ptr_mc_1 < ptr_mc_0 + cz; ptr_mb_1++, ptr_mc_1++)
            {
                *ptr_mc_1 += buf_ma * *ptr_mb_1;
            }

            ptr_ma_1++;
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulacc_test_v1() PASSED\n");
#endif

    return ierr;
}


int stblas_3e_mxmxmulacc_test_v2(const int rdr, const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double *ma, const size_t lda, const double *mb, const size_t ldb, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi, zj, zk;

    double buf_mc;

#ifdef TRACE_INFO
    printf("[F] stblas_3e_mxmxmulacc_test_v2()\n");
#endif

    for (zj = 0; zj < cz; zj++)
    {
        for (zi = 0; zi < rz; zi++)
        {
            buf_mc = 0;

            for (zk = 0; zk < kz; zk++)
            {
                buf_mc += ma[zk + zi * lda] * mb[zj + zk * ldb];
            }

            mc[zj + zi * ldc] = buf_mc;
        }
    }

#ifdef TRACE_INFO
    printf("    stblas_3e_mxmxmulacc_test_v2() PASSED\n");
#endif

    return ierr;
}


int stblas_gemm_test(const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double alpha, const double *ma, const size_t lda, const double *mb, const size_t ldb, const double beta, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi, zj, zk;

    size_t a_rz, a_cz, b_rz;

    int flg_ta, flg_tb;
    double temp;

#ifdef TRACE_INFO
    printf("[F] stblas_gemm_test()\n");
#endif

    flg_ta = trna;
    flg_tb = trnb;

    if (trna == MTBLAS_TRANS_NONE)
    {
        a_rz = rz;
        a_cz = kz;
    }
    else
    {
        a_rz = kz;
        a_cz = rz;
    }

    if (trnb == MTBLAS_TRANS_NONE)
    {
        b_rz = kz;
    }
    else
    {
        b_rz = cz;
    }

    if (trna != MTBLAS_TRANS_NONE && trna != MTBLAS_TRANS && trna != MTBLAS_TRANS_CONJ)
    {
        ierr = ERRNO_UNKNOWN_TYPE;
    }
    else if (trnb != MTBLAS_TRANS_NONE && trnb != MTBLAS_TRANS && trnb != MTBLAS_TRANS_CONJ)
    {
        ierr = ERRNO_UNKNOWN_TYPE;
    }
    else if (rz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (cz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (kz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (lda < 1 || lda < a_rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }
    else if (ldb < 1 || ldb < b_rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }
    else if (ldc < 1 || ldc < rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }

    if (ierr != ERRNO_UNDEFINED)
    {
        printf("ERROR!\n");
        goto EXIT_STBLAS_GEMM_TEST;
    }

    if ((rz == 0 || cz == 0) || (alpha == 0 || kz == 0) && (beta == 1.0))
    {
        goto EXIT_STBLAS_GEMM_TEST;
    }

    if (alpha == 0)
    {
        if (beta == 0)
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    mc[zi + zj * ldc] = 0;
                }
            }
        }
        else
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    mc[zi + zj * ldc] = beta * mc[zi + zj * ldc];
                }
            }
        }

        goto EXIT_STBLAS_GEMM_TEST;
    }

    if (flg_tb == MTBLAS_TRANS_NONE)
    {
        if (flg_ta == MTBLAS_TRANS_NONE)
        {
            for (zj = 0; zj < cz; zj++)
            {
                if (beta == 0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = 0;
                    }
                }
                else if (beta != 1.0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = beta * mc[zi + zj * ldc];
                    }
                }
                for (zk = 0; zk < kz; zk++)
                {
                    if (mb[zk + zj * ldb] != 0)
                    {
                        temp = alpha * mb[zk + zj * ldb];
                        for (zi = 0; zi < rz; zi++)
                        {
                            mc[zi + zj * ldc] = mc[zi + zj * ldc] + temp * ma[zi + zk * lda];
                        }
                    }
                }
            }
        }
        else
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    temp = 0;
                    for (zk = 0; zk < kz; zk++)
                    {
                        temp = temp + ma[zk + zi * lda] * mb[zk + zj * ldb];
                    }
                    if (beta == 0)
                    {
                        mc[zi + zj * ldc] = alpha * temp;
                    }
                    else
                    {
                        mc[zi + zj * ldc] = alpha * temp + beta * mc[zi + zj * ldc];
                    }
                }

            }
        }
    }
    else
    {
        if (flg_ta == MTBLAS_TRANS_NONE)
        {
            for (zj = 0; zj < cz; zj++)
            {
                if (beta == 0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = 0;
                    }
                }
                else
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = beta * mc[zi + zj * ldc];
                    }
                }
                for (zk = 0; zk < kz; zk++)
                {
                    if (mb[zj + zk * ldb] != 0)
                    {
                        temp = alpha * mb[zj + zk * ldb];
                        for (zi = 0; zi < rz; zi++)
                        {
                            mc[zi + zj * ldc] = mc[zi + zj * ldc] + temp * ma[zi + zk * lda];
                        }
                    }
                }
            }
        }
        else
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    temp = 0;

                    for (zk = 0; zk < kz; zk++)
                    {
                        temp += ma[zk + zi * lda] * mb[zj + zk * ldb];
                    }

                    if (beta == 0)
                    {
                        mc[zi + zj * ldc] = alpha * temp;
                    }
                    else
                    {
                        mc[zi + zj * ldc] = alpha * temp + beta * mc[zi + zj * ldc];
                    }
                }
            }
        }
    }

    EXIT_STBLAS_GEMM_TEST:
#ifdef TRACE_INFO
    printf("    stblas_gemm_test() PASSED\n");
#endif

    return ierr;
}


int stblas_gemm_test_v1(const int trna, const int trnb, const size_t rz, const size_t cz, const size_t kz, const double alpha, const double *ma, const size_t lda, const double *mb, const size_t ldb, const double beta, double *mc, const size_t ldc)
{
    int ierr = ERRNO_UNDEFINED;
    size_t zi, zj, zk;

    size_t a_rz, a_cz, b_rz;

    int flg_ta, flg_tb;
    double temp;

#ifdef TRACE_INFO
    printf("[F] stblas_gemm_test_v1()\n");
#endif

    flg_ta = trna;
    flg_tb = trnb;

    if (trna == MTBLAS_TRANS_NONE)
    {
        a_rz = rz;
        a_cz = kz;
    }
    else
    {
        a_rz = kz;
        a_cz = rz;
    }

    if (trnb == MTBLAS_TRANS_NONE)
    {
        b_rz = kz;
    }
    else
    {
        b_rz = cz;
    }

    if (trna != MTBLAS_TRANS_NONE && trna != MTBLAS_TRANS && trna != MTBLAS_TRANS_CONJ)
    {
        ierr = ERRNO_UNKNOWN_TYPE;
    }
    else if (trnb != MTBLAS_TRANS_NONE && trnb != MTBLAS_TRANS && trnb != MTBLAS_TRANS_CONJ)
    {
        ierr = ERRNO_UNKNOWN_TYPE;
    }
    else if (rz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (cz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (kz < 0)
    {
        ierr = ERRNO_UNDEFINED;
    }
    else if (lda < 1 || lda < a_rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }
    else if (ldb < 1 || ldb < b_rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }
    else if (ldc < 1 || ldc < rz)
    {
        ierr = ERRNO_OUT_OF_RANGE_SIZE;
    }

    if (ierr != ERRNO_UNDEFINED)
    {
        printf("ERROR!\n");
        goto EXIT_STBLAS_GEMM_TEST;
    }

    if ((rz == 0 || cz == 0) || (alpha == 0 || kz == 0) && (beta == 1.0))
    {
        goto EXIT_STBLAS_GEMM_TEST;
    }

    if (alpha == 0)
    {
        if (beta == 0)
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    mc[zi + zj * ldc] = 0;
                }
            }
        }
        else
        {
            for (zj = 0; zj < cz; zj++)
            {
                for (zi = 0; zi < rz; zi++)
                {
                    mc[zi + zj * ldc] *= beta;
                }
            }
        }

        goto EXIT_STBLAS_GEMM_TEST;
    }

    if (flg_tb == MTBLAS_TRANS_NONE)
    {
        if (flg_ta == MTBLAS_TRANS_NONE)
        {
            for (zj = 0; zj < cz; zj++)
            {
                if (beta == 0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = 0;
                    }
                }
                else if (beta != 1.0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = beta * mc[zi + zj * ldc];
                    }
                }
                for (zk = 0; zk < kz; zk++)
                {
                    if (mb[zk + zj * ldb] != 0)
                    {
                        temp = alpha * mb[zk + zj * ldb];
                        for (zi = 0; zi < rz; zi++)
                        {
                            mc[zi + zj * ldc] = mc[zi + zj * ldc] + temp * ma[zi + zk * lda];
                        }
                    }
                }
            }

            /*
            if (beta == 0)
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = 0;
                    }

                    for (zk = 0; zk < kz; zk++)
                    {
                        if (mb[zk + zj * ldb] != 0)
                        {
                            temp = alpha * mb[zk + zj * ldb];
                            for (zi = 0; zi < rz; zi++)
                            {
                                mc[zi + zj * ldc] += temp * ma[zi + zk * lda];
                            }
                        }
                    }
                }
            }
            else if (beta != 1.0)
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] *= beta;
                    }

                    for (zk = 0; zk < kz; zk++)
                    {
                        if (mb[zk + zj * ldb] != 0)
                        {
                            temp = alpha * mb[zk + zj * ldb];
                            for (zi = 0; zi < rz; zi++)
                            {
                                mc[zi + zj * ldc] += temp * ma[zi + zk * lda];
                            }
                        }
                    }
                }
            }
            else
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zk = 0; zk < kz; zk++)
                    {
                        if (mb[zk + zj * ldb] != 0)
                        {
                            temp = alpha * mb[zk + zj * ldb];
                            for (zi = 0; zi < rz; zi++)
                            {
                                mc[zi + zj * ldc] += temp * ma[zi + zk * lda];
                            }
                        }
                    }
                }
            }
            */
        }
        else
        {
            if (beta == 0)
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        temp = 0;
                        for (zk = 0; zk < kz; zk++)
                        {
                            temp += ma[zk + zi * lda] * mb[zk + zj * ldb];
                        }
                        mc[zi + zj * ldc] = alpha * temp;
                    }
                }
            }
            else if (beta == 1.0)
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        temp = 0;
                        for (zk = 0; zk < kz; zk++)
                        {
                            temp += ma[zk + zi * lda] * mb[zk + zj * ldb];
                        }
                        mc[zi + zj * ldc] += alpha * temp;
                    }
                }
            }
            else
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        temp = 0;
                        for (zk = 0; zk < kz; zk++)
                        {
                            temp += ma[zk + zi * lda] * mb[zk + zj * ldb];
                        }
                        mc[zi + zj * ldc] *= beta;
                        mc[zi + zj * ldc] += alpha * temp;
                    }
                }
            }
        }
    }
    else
    {
        if (flg_ta == MTBLAS_TRANS_NONE)
        {
            for (zj = 0; zj < cz; zj++)
            {
                if (beta == 0)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] = 0;
                    }
                }
                else
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        mc[zi + zj * ldc] *= beta;
                    }
                }
                for (zk = 0; zk < kz; zk++)
                {
                    if (mb[zj + zk * ldb] != 0)
                    {
                        temp = alpha * mb[zj + zk * ldb];
                        for (zi = 0; zi < rz; zi++)
                        {
                            mc[zi + zj * ldc] += temp * ma[zi + zk * lda];
                        }
                    }
                }
            }
        }
        else
        {
            if (beta == 0)
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        temp = 0;
                        for (zk = 0; zk < kz; zk++)
                        {
                            temp += ma[zk + zi * lda] * mb[zj + zk * ldb];
                        }

                        mc[zi + zj * ldc] = alpha * temp;
                    }
                }
            }
            else
            {
                for (zj = 0; zj < cz; zj++)
                {
                    for (zi = 0; zi < rz; zi++)
                    {
                        temp = 0;
                        for (zk = 0; zk < kz; zk++)
                        {
                            temp += ma[zk + zi * lda] * mb[zj + zk * ldb];
                        }
                        mc[zi + zj * ldc] *= beta;
                        mc[zi + zj * ldc] += alpha * temp;
                    }
                }
            }
        }
    }

    EXIT_STBLAS_GEMM_TEST:
#ifdef TRACE_INFO
    printf("    stblas_gemm_test_v1() PASSED\n");
#endif

    return ierr;
}


int test_3e_mtxmtxmul(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    mat_df_t ma;
    mat_df_t mb;
    mat_df_t mc;
    mat_df_t md;

    struct timespec tp_b;  /* Time of test begin */
    struct timespec tp_e;  /* Time of test end */

    ierr = mat_df_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_df_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);
    ierr = mat_df_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    printf("PROPER_R_SIZE_RUN: %d, PROPER_C_SIZE_RUN: %d\n", PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN);
    printf("Test: stblas_2e_mxselfscal()\n");
    printf("Test: stblas_3e_mxmxmulacc()\n");
/*
    printf("Row-Major:\n");
    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 1; ii++)
    //for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 1; ij++)
        //for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            ierr = stblas_2e_mxselfscal(MTBLAS_ROW_MAJOR, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, 0, md.data, md.row_size);

            //ierr = stblas_3e_mxmxmulacc(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, ma.data, ma.row_size, mb.data, mb.row_size, md.data, md.row_size);
            //ierr = stblas_3e_mxmxmulacc_test(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, ma.data, ma.row_size, mb.data, mb.row_size, md.data, md.row_size);
            ierr = stblas_3e_mxmxmulacc_test_v1(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, ma.data, ma.row_size, mb.data, mb.row_size, md.data, md.row_size);
            //ierr = stblas_3e_mxmxmulacc_test_v2(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, ma.data, ma.row_size, mb.data, mb.row_size, md.data, md.row_size);

        #ifdef DEBUG_INFO
            printf("%d: stblas_3e_mxmxmulacc()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (stblas_2e_mxselfscal() + stblas_3e_mxmxmulacc())\n", clock_difftime(&tp_b, &tp_e));

    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 1; ii++)
    //for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 1; ij++)
        //for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            cblas_dgemm(CblasRowMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, 1.0, ma.data, ma.row_size, mb.data, mb.row_size, 0, md.data, md.row_size);

        #ifdef DEBUG_INFO
            printf("%d: cblas_dgemm()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (cblas_dgemm())\n", clock_difftime(&tp_b, &tp_e));
*/
    printf("Column-Major:\n");
/*
    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 2; ii++)
    {
        for (ij = 0; ij < 2; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            ierr = stblas_2e_mxselfscal(MTBLAS_COL_MAJOR, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, 0, md.data, md.row_size);

            ierr = stblas_3e_mxmxmulacc(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, ma.data, ma.col_size, mb.data, mb.col_size, md.data, md.col_size);

        #ifdef DEBUG_INFO
            printf("%d: stblas_3e_mxmxmulacc()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (stblas_2e_mxselfscal() + stblas_3e_mxmxmulacc())\n", clock_difftime(&tp_b, &tp_e));
*/
    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 2; ii++)
    {
        for (ij = 0; ij < 2; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            ierr = stblas_gemm_test(flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, 1.0, ma.data, ma.col_size, mb.data, mb.col_size, 0, md.data, md.col_size);

        #ifdef DEBUG_INFO
            printf("%d: stblas_gemm_test()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (stblas_gemm_test())\n", clock_difftime(&tp_b, &tp_e));

    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 2; ii++)
    {
        for (ij = 0; ij < 2; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            ierr = stblas_gemm_test_v1(flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, 1.0, ma.data, ma.col_size, mb.data, mb.col_size, 0, md.data, md.col_size);

        #ifdef DEBUG_INFO
            printf("%d: stblas_gemm_test_v1()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (stblas_gemm_test_v1())\n", clock_difftime(&tp_b, &tp_e));

    clock_gettime(CLOCK_REALTIME, &tp_b);
    cnt = 1;
    for (ii = 0; ii < 2; ii++)
    {
        for (ij = 0; ij < 2; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);

            cblas_dgemm(CblasColMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, 1.0, ma.data, ma.col_size, mb.data, mb.col_size, 0, md.data, md.col_size);

        #ifdef DEBUG_INFO
            printf("%d: cblas_dgemm()\n", cnt);
            mat_df_stdout(&md);
        #endif
            cnt++;
        }
    }
    clock_gettime(CLOCK_REALTIME, &tp_e);
    printf("Time: %10.3lf (cblas_dgemm())\n", clock_difftime(&tp_b, &tp_e));

    ierr = mat_df_delete(&ma);
    ierr = mat_df_delete(&mb);
    ierr = mat_df_delete(&mc);
    ierr = mat_df_delete(&md);

    return ierr;
}
